# x0s = [[1.2; 1.2], [-1.2; 1]];
# for i = 1:size(x0s, 2)
# x0 = x0s(:, i);
# [xk, fk, gk, Bk] = BFGS(@rosenbrock, @rosenbrock_gradient, @rosenbrock_hessian, x0);
# fprintf('Initial point: (%f, %f)\n', x0(1), x0(2));
# fprintf('Solution point: (%f, %f)\n', xk(1), xk(2));
# fprintf('Function value: %f\n', fk);
# fprintf('Gradient norm: %f\n', norm(gk));
# fprintf('Frobenius norm of Bk - Hessian: %f\n\n', norm(Bk - rosenbrock_hessian(xk)));
# end
#
#
# function [xk, fk, gk, Bk] = BFGS(f, g, H, x0)
# max_iter = 10000;
# epsilon = 1e-6;
# alpha = 0.5;
# beta = 0.5;
# xk = x0;
# gk = g(xk);
# Bk = eye(length(x0));
# k = 0;
# while norm(gk) > epsilon && k < max_iter
#     pk = -Bk * gk;
#     rho = alpha;
#     while f(xk + rho * pk) > f(xk) + alpha * rho * gk' * pk
#     rho = beta * rho;
# end
# sk = rho * pk;
# yk = g(xk + sk) - gk;
# if sk' * yk > 0
# Bk = (eye(length(x0)) - sk * yk' / (yk' * sk)) * Bk * (eye(length(x0)) - yk * sk' / (yk' * sk)) + sk * sk' / (yk' * sk);
# end
# xk = xk + sk;
# gk = g(xk);
# k = k + 1;
# end
# fk = f(xk);
# end
#
# function val = rosenbrock(x)
# val = 100*(x(2)-x(1)^2)^2 + (1-x(1))^2;
# end
#
# function grad = rosenbrock_gradient(x)
# grad = [-400*x(1)*(x(2)-x(1)^2)-2*(1-x(1)); 200*(x(2)-x(1)^2)];
# end
#
# function hess = rosenbrock_hessian(x)
# hess = [1200*x(1)^2-400*x(2)+2, -400*x(1); -400*x(1), 200];
# end